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ABSTRACT 

The transiting exoplanet HD 149026b is an important case for theories of planet formation and planetary 
structure, because the planet's relatively small size has been interpreted as evidence for a highly metal-enriched 
\ composition. We present observations of 4 transits with the Near Infrared Camera and Multi-Object Spectrom- 

eter on the Hubble Space Telescope within a wavelength range of 1.1-2.0 /im. Analysis of the light curve 
gives the most precise estimate yet of the stellar mean density, p+ = 0.497±[j g cm" 3 . By requiring agreement 
CN) ! between the observed stellar properties (including p+) and stellar evolutionary models, we refine the estimate 

■ of the stellar radius: = 1. 541^ ^2 ^o- We also find a deeper transit than has been measured at optical and 
q_) I mid-infrared wavelengths. Taken together, these findings imply a planetary radius of R p = 0.813^'o|g /?j up , 

■ which is larger than earlier estimates. Models of the planetary interior still require a metal-enriched composi- 
tion, although the required degree of metal enrichment is reduced. It is also possible that the deeper NICMOS 
transit is caused by wavelength-dependent absorption by constituents in the planet's atmosphere, although sim- 

,—, . pie model atmospheres do not predict this effect to be strong enough to account for the discrepancy. We use 

PLh ' the 4 newly-measured transit times to compute a refined transit ephemeris. 

P_] . Subject headings: stars: individual (HD 149026) — techniques: photometric — stars: planetary systems — 

(-h stars: fundamental parameters 

Pk 

O 1 1- INTRODUCTION 

■ Since its discovery by Sato et al. (2005), HD 149026b has been one of the most closely scrutinized planets outside the Solar 
system. It is a close-in gas giant, orbiting a G star with a period of only 2.5 d. Observations of transits (Sato et al. 2005, 

1— 1 , Charbonneau et al. 2006, Winn et al. 2008b, Nutzman et al. 2008), in combination with observations of radial-velocity variations 
of the parent star (Sato et al. 2005), have shown that the planet has approximately Saturn's mass but is considerably denser, 
despite the intense irradiation from the parent star that should inflate the planet and lower its density. There is consensus among 
theorists that the reason for the "shrunken radius" is a highly metal-enriched composition, although the total metal mass, its 
distribution within the planet, and the reason for the enrichment are debated [Sato et al. (2005), Fortney et al. (2006), Ikoma et 
al. (2006), Broeg & Wuchterl (2007), Burrows et al. (2007)]. The total metal mass, for example, ranges from 60 M e to 1 14 M® 
among the possible models. The latter estimate would represent 80% of the total mass of the planet. 

The planet's outer atmosphere is also of interest, given the possibly unusual composition and the strong heating from the parent 
star. Models by Fortney et al. (2005) indicated the possibility of a very hot stratosphere as a result of gaseous TiO and VO 
opacity. By using the Spitzer Space Telescope to observe a planetary occultation, Harrington et al. (2007) found the planet's 
8 /im brightness temperature to be much larger than the temperature that one would expect based on thermal equilibrium with the 
incident stellar radiation. This may be the result of the predicted TiO and VO heating, although the details of whether and where 
these absorbers actually condense in the atmosphere are not yet understood (Fortney et al. 2005, Harrington et al. 2007). 

Fundamental to all these discussions are the measurements of the mass and radius of HD 149026b. These measurements are 
limited by the uncertainties in the stellar mass and radius. One way to improve the situation is to observe transits with greater 
photometric precision than has been possible before. As shown by Seager & Mallen-Ornelas (2003), with a good light curve and 
Kepler's third law, one may determine the stellar mean density. If the mean density is known precisely enough, it is a key con- 
straint that can be combined with the other stellar observables (parallax, apparent magnitude, effective temperature, metallicity, 
etc.) and stellar-evolutionary models to determine the stellar mass and radius. This technique has been put into practice for many 
other systems [see, e.g., Sozzetti et al. (2007), Holman et al. (2007), Torres et al. (2008)] but never to advantage for HD 149026b 
because of the limited precision of prior determinations of p+ (Winn et al. 2008b, Nutzman et al. 2008). Observers must cope with 
the small transit depth of 2.5 mmag (smaller than any other transiting exoplanet by a factor of two) and the paucity of suitable 
comparison stars within the field of view of most telescopes. 

In this paper we present observations of transits of HD 149026b with the Near Infrared Camera and Multi-Object Spectrometer 
[NICMOS, Thompson (1992)] on board the Hubble Space Telescope (HST). We chose this instrument because a high precision 
in relative photometry is possible even without using comparison stars (Gilliland 2006, Swain et al. 2008) and because the 
reduced stellar limb-darkening at near-infrared wavelengths is advantageous for the light-curve analysis (Carter et al. 2008, Pal 
et al. 2008). We have organized this paper as follows. In §|2]we describe the observations and data reductions leading to the 
final photometric time series. In §|3]we describe our photometric model and the results of the NICMOS light-curve analysis. In 
§ |U we describe how the light-curve results were incorporated into stellar-evolutionary models to determine the parameters of 
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the HD 149026 system. In § [5] the light analysis is repeated using not only the NICMOS data but also the most precise light 
curves that have been obtained at optical and mid-infrared wavelengths. In § [6] all the available transit times are analyzed to 
produce a refined transit ephemeris and to search for possible period variations that could be indicative of additional bodies in the 
HD 149026 system (Holman & Murray 2005, Agol et al. 2005, Ford & Holman 2007). Finally, in §0 we discuss the possible 
implications of our observations and analysis. 

2. OBSERVATIONS AND REDUCTIONS 

We observed HD 149026 on four occasions ("visits" in HST parlance) when transits were predicted to occur, on 2007 Dec 22, 
2007 Dec 24, 2008 Feb 08, and 2008 Mar 20. Each visit consisted of five orbits spanning a transit. Between each pair of orbits is 
an observing gap of approximately 45 minutes, when HST's vision is blocked by the Earth. The visits were scheduled in such a 
manner that the combined data set provides complete phase coverage of the transit, including redundant coverage of the critical 
ingress and egress phases. In particular, visits 1 and 3 covered the ingress phase, and visit 2 covered both ingress and egress 
phases. Visit 4 captured the beginning of egress. 

We used Camera 3 of the NICMOS detector, a 256 x 256 HgCdTe array with a field of view of 51.2." x 51.2.". We used the 
G141 grism filter, which is centered at 1 A/im, spans 0.8 [im and is roughly equivalent to H band. An exposure was obtained every 
13 s. The camera was operated in "MULTIACCUM" mode, wherein five nondestructive readouts are recorded during a single 
exposure, and the first readout is subtracted from the final readout. After accounting for overheads, the effective integration time 
was 4 s per exposure. We deliberately defocused the instrument to give a full- width at half-maximum (FWHM) of approximately 
5 pixels in the cross-dispersion direction. This was done for two reasons: firstly, when focused, camera 3 undersamples the 
point-spread-function (PSF) of point sources; and secondly, the detector pixels exhibit intra-pixel sensitivity variations as large 
as 30%. Defocusing the images causes the PSF to be well-sampled and averages over the intra-pixel sensitivity variations. 

Approximately 220 exposures were collected during each HST orbit. Experience with HST has shown that photometric 
stability is relatively poor during the first orbit of a given visit. Our observations were scheduled under the assumption that the 
first orbit from each visit would not be utilized, and indeed we ended up omitting the first-orbit data from our analysis. At the 
start of each visit, we obtained a single non-dispersed image using a narrow filter centered at 1 .66 fim in order to establish the 
pixel position corresponding to zero dispersion. We then adopted the relation from the HST Data Handbook for NICMOS 4 , 

X(Ax) = -0.007992 Ax+ 1.401 (1) 

where X(Ax) is the wavelength (measured in /im), and Ax is the x coordinate (measured in pixels) relative to the center of the 
undispersed image. 

For completeness, we performed the steps of flat-fielding, background subtraction, and pixel flagging, as described below; 
however, it is noteworthy that these steps in the data reduction made very little difference in the aperture photometry or in the 
final results. Flat-field correction for grism images is not straightforward and is not done as part of the standard NICMOS 
pipeline reductions, because the appropriate flat field depends both upon wavelength and upon the position of the source in 
the non-dispersed image. To accomplish the flat-field correction, we obtained seven flat fields, each using a narrow bandpass 
within the G141 bandpass, and fitted the data for each pixel with a quadratic function of wavelength C[X;x,y]. We then applied 
C[X(Ax);x,y] as a multiplicative correction to each of our science images, using the wavelength-coordinate relation A(Ajt) that 
was determined from the single non-dispersed image. The background level was estimated in each image based on the counts in a 
relatively clean region of the detector (away from the spectral trace) and subtracted from the entire image. To identify bad pixels, 
all images from a given orbit were used to create a time series of counts specific to each pixel. Pixels showing an anomalously 
large variance were flagged. The list of flagged pixels was appended to the list of hot or cold pixels that were identified in 
the standard NICMOS pipeline reductions, and the values of all of those bad pixels were replaced by interpolated values of the 
neighboring good pixels. 

Aperture photometry was performed on the first-order spectrum, using a simple sum of the counts within a rectangular box 
centered on the spectral trace. The box had a width of 20 pixels in the cross-dispersion direction (the y direction), which was four 
times the FWHM of the PSF. The box had a length of 120 pixels in the dispersion direction (the x direction), which was long 
enough the capture the entire first-order spectrum. 

At this stage the data had been reduced to a single number per image: the total number of counts in the aperture (the "flux"). 
We examined the resulting time series. As expected, the data collected during the first orbits of each visit showed flux variations 
that were both larger in amplitude and different in their time-dependence than the variations observed in subsequent orbits. The 
first-orbit data were excluded from subsequent analysis. In addition, we excluded the 10 exposures near the beginning and end of 
each orbit sequence, because they showed strong flux variations that are probably due to "Earth shine." After these exclusions, 
there remained 800, 820, 800, and 792 good data points in visits 1, 2, 3, and 4, respectively. 

3. NICMOS LIGHT-CURVE ANALYSIS 

Fig.Q]shows the time series of the aperture-summed flux, after dividing by the mean flux. In each panel, the zero point of the 
x-axis is the expected mid-transit time. The flux decrement during the transit is identifiable, but this decrement is superimposed 
on at least two other sources of variability: orbit-to-orbit discontinuities, and smooth intra-orbital variability showing a consistent 
pattern among all orbits of a given visit. The intra-orbital variability has been seen by all other investigators attempting precise 
HST photometry of single bright stars, since the pioneering work by Brown et al. (2001), and the flux discontinuities have been 
seen by other investigators using NICMOS (see, e.g., Swain et al. 2008). The origins of these systematic effects have not been 
established. The orbit-to-orbit consistency of the smooth variations suggests a phenomenon that is a function of the phase of 
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the telescope's orbit around the Earth, such as thermal cycling or scattered light. The discontinuities between orbits suggest a 
non-repeating event associated with the re-aquisition of the target star after each Earth occultation, such as pointing changes or 
positional shifts of the grism filter. 

Ideally, the underlying physical processes giving rise to these systematic effects could be ascertained, and this understanding 
would lead to either the recognition of an improved method for deriving the photometric signal or a physical model that could be 
used to correct the aperture-summed flux. Given that we do not yet have such knowledge, what can be done? The intra-orbital 
variations are very well-described by a smooth function of the HST orbital phase; following other investigators we used a smooth 
function with several adjustable parameters as an ad hoc model for this variation. The parameters of this model are fairly well 
constrained by the out-of-transit data, for which all variations are assumed to be systematic effects. 

The inter-orbital discontinuities are more problematic. To investigate the systematic effects, we examined the spectral trace on 
each image. Specifically we computed the flux-weighted mean y position as a function of x, giving a curve y(x) representing the 
centroid of the spectral trace (the "footprint" of the spectrum on the detector). We also estimated the orientation of the spectral 
trace relative to the detector edges, by performing a linear fit to the previously calculated function y(x). figure |2] shows the 
results. Within a single orbit, the position and orientation of the spectral trace are relatively constant, as compared to the larger 
movements that are observed between orbits. The largest variations of the spectral trace (its position and width) seem to coincide 
with the largest discontinuities in the flux time series. 

Given the correlations that are observed between the properties of the spectral trace and the aperture-summed flux, the ap- 
proach taken by Swain et al. (2008) and other investigators is to "decorrelate" the flux against a number of measured parameters 
("state variables") such as the mean y position, cross-dispersion width, orientation angle, and so forth. One way to achieve this 
decorrelation is to fit linear functions of the state variables to the out-of-transit data, for which all time variations are expected to 
be due to the systematic effects. Then the best-fitting parameters are used to correct all of the data. Alternatively, one could fit 
for the linear functions of the state variables simultaneously with the parameters describing the transit light curve. 

We attempted both of these procedures and found that while they do reduce the amplitude of the systematic effects, they 
still leave highly significant systematic variations. We also find this procedure to be undesirable because it is not clear which 
parameters to include in the fit; because the fitted parameters are highly correlated (the state variables do not vary independently); 
and because we have no justification for the assumption of a linear function for any of these parameters, without an understanding 
of the underlying physical effect. An example of a possibly relevant physical effect that would not necessarily be described by 
a linear function is intra-pixel sensitivity variation, which could lead to a function that is periodic in the pixel coordinates of the 
spectral trace. 

We attempted to fit numerous physically-motivated models (based on the premise of intra-pixel sensitivity variations, among 
others), and did not find any such model that provided a good fit to the data while also having only a few, nondegenerate adjustable 
parameters. Ultimately we gave up on attempting to correct the intra-orbital discontinuities based on a priori information. Instead 
we included in our model an adjustable multiplicative factor specific to each orbit. This might seem devastating to the goal of 
analyzing the transit light curve, but this is not so. It means that we cannot make use of the relative flux between orbits to derive 
the transit depth; but as we will show, the truly precious information is in the duration of ingress or egress, which is much less 
vulnerable to the problem of flux discontinuities. In addition, four of the orbits spanned a full ingress or egress. Data from those 
four orbits does provide useful information about the transit depth, because the discontinuities appear between orbits and not 
within orbits. 

All together, our model for the flux variation due to systematic effects is 

/sys(O = /o'{l+^(O + C , 1 [0(f)] 2 + ^[0(f)] 3 } , (2) 

where the v index specifies the visit number (1^4), the o index specifies the orbit number (1^4) within each visit after omitting the 
first orbit, the 16 numbers f£ are the multiplicative factors describing the inter-orbital discontinuities, <fi(t) is the HST orbital phase 
at time t, and the 12 numbers c] (3 per visit) are constants specifying a polynomial function of tj> that describes the intra-orbital 
variation. The HST orbital phase was defined as 

(t-(t)) mod Phst 

m ^^2L ™I (3) 

"est 

where Phst = 1 -5975 hours is the orbital period of the HST around the Earth and (t) is the midpoint of each orbit's observations. 
The choice of a polynomial, as opposed to some other smoothly varying function, was arbitrary. We also tried using sinusoidal 
functions with an angular frequency of 2-k/Phst, with no significant differences in any of the results described below. 

For the transit model, we used the analytic formulas of Mandel & Agol (2002). Our parameters were the planet-to-star radius 
ratio (R p /R+), the cosine of the orbital inclination (cos/), the semimajor axis in units of the stellar radius (a/R+), the mid-transit 
time t*, and the two coefficients u\ and ui of a quadratic limb-darkening law, 

^ = l- Ml (l- A1 )- M2 (l-/i) 2 (4) 

where /i is the cosine of the angle between the observer and the normal to the stellar surface and 1^ is the specific intensity as a 
function of /i. We allowed ui and U2 to vary freely, subject to the conditions u\ + U2 < 1, u\ + U2 > 0, and u\ > 0. These conditions 
require the brightness profile to be everywhere positive and monotonically decreasing from limb to center. In practice, the fitting 
parameters were actually 

u[ = mi cos 40°- ii2 sin 40° (5) 
u' 2 = ui sin 40° + M 2 cos 40° (6) 
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FIG. 1 . — NICMOS photometry (1.1-2.0 fim) of HD 149026b of 4 transits, with interruptions due to Earth occultations. Plotted are the results of simple aperture 
photometry. The observed variations are a combination of the transit signal and systematic effects (intra-orbital variations and inter-orbital discontinuities). The 
solid curve is the best-fitting model that accounts for both the transit signal and systematic effects. 

because u\ and u' 2 are weakly correlated, unlike U\ and m_ (Pal 2008). In computing the transit light curve we assumed the orbit to 
be circular, consistent with the findings of Sato et al. (2005) and Madhusudhan & Winn (2008). We held the orbital period fixed 
at the value P = 2.87588 days based on the results of Winn et al. (2008b). Here the period is used only to relate the measured 
transit durations and a/R+. The fractional error in P is approximately 10 4 times smaller than the fractional error in a/R* and is 
safely ignored (Carter et al. 2008). 

Our complete model of the photometric time series was the product of the transit model and f sys (t) (Eq. [2j. We fitted simulta- 
neously for the parameter set R p /R+, cos;, a/R+, u\, u' 2 , {t^}, {/„}, and {cj'}). The polynomial describing intra-orbital variations 
was specific to each visit, and the flux discontinuities were specific to each orbit, but the transit parameters were required to be 
consistent across all orbits and visits. We performed a least-squares fit to the unbinned data using a box-constrained Levenberg- 
Marquardt algorithm (Levenberg 1944, Marquardt 1963, Lourakis 2004) utilizing the Jacobian calculation of Pal (2008). Box 
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FIG. 2. — Illustration of inter-orbital variations of the spectral trace. The solid curves are the flux-weighted mean y position of the first-order spectrum as a 
function of x. Overplotted are linear fits to v(.v). The inset figure shows the measured light curve after dividing out the intra-orbital variations correlated with HST 
orbital phase. The largest rotation of the spectral trace (at the fifth orbit) coincides with the largest discontinuity in the flux time series. 

constraints were needed to enforce the restrictions on the limb darkening parameters u\ and ui. The goodness-of-fit statistic was 



X 



=1 i=l 



(7) 



where P mo( fS) the calculated flux at the time of the i th data point during visit v, fl bs (f) is the i th flux measurement during visit 
v, N v is the number of data points in visit v, and a v was assumed to be a constant at this step. The solid curve in Fig. Q] shows 
the best-fitting model. The root-mean-square (rms) residual between the data and the best-fitting model is 440 parts per million 
(ppm). This is approximately 2.2 times the expected noise level calculated within the NICMOS calibration pipeline (which is 
dominated by photon noise). Fig.[3]shows histograms of the residuals for all the data and for each visit individually. The residuals 
are not Gaussian; the flattened peak in the histograms indicates a nonzero kurtosis 5 . 

Evidently the noise is not photon-limited, and is not Gaussian, but at least it does not appear to be strongly correlated in time. 
We assessed the degree of temporal correlations ("red noise") in two ways. First, we binned the residuals in time by a factor yV 
ranging from 1 to 100, and calculated the standard deviation ijn of the binned residuals. The results are shown in Fig. [4] They 
follow closely the expectation of independent random numbers, ojv = &i N~ x l 2 [M/(M- l)] 1 / 2 , where M is the number of bins. 
Second, we calculated the Allan (1964) variance a\{l) of the residuals, defined as 



1 



2(N+\-2l) 



N-21 I . l-\ 
i=0 V 7=0 



i+j r i+ j+i 



(8) 



where denotes the residual of the kih data point, is the number of data points, and / is the lag. The Allan variance is 
commonly used in the time metrology literature to assess 1 // noise. For independent residuals, one expects a\(J) ~ a\(Q))/l. The 
results for our data, also shown in Fig.@] satisfy this expectation. There is no readily identifiable time-correlated component in 
the time series of residuals. 



One may wonder about the effect of the apparent non-Gaussianity of the noise, shown in Fig. 3. To investigate this issue we used an Edgeworth expansion to 
create a new x 2 -like statistic that accounts for the skewness and kurtosis of the residuals (see, e.g., Amendola et al. 1996). Using this different fitting statistic, we 
found that the best-fit parameter values were unchanged. This was not surprising because bias is expected to arise from skewness (as opposed to kurtosis) and the 
skewness of the residuals is very small. However, the confidence intervals are affected by the kurtosis. We found that accounting for kurtosis leads to error bars 
that are smaller than the error bars quoted here, but only by a modest amount (< 20%). For simplicity, the results quoted in this paper are based on the standard 
X 2 statistic given in Eq. (7). 
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FIG. 3. — Histograms of the residuals between the data and the best-fitting model. The solid line is the histogram based on all the data. The other lines are for 
data specific to visit 1 (dotted), visit 2 (dashed), visit 3 (dash-dot), and visit 4 (dash-dot-dot-dot). 
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FIG. 4. — Assessment of correlated noise. Left panel: The rms of time-binned residuals, as a function of bin size, for each of the 4 visits. Right panel: The 
Allan variance of the residuals, as a function of lag, for each of the 16 orbits. The dotted lines are the results of the calculations based on the data, and the solid 
lines show the expected trend if the noise were uncorrelated. 

Plotted in Fig.[5]is the measured flux after dividing by the optimized function f sys (t). This represents our best effort to correct 
for the systematic effects. In Fig. [6] we show the results of combining the data from all the visits (after correcting for systematic 
effects) into a single transit light curve. In this composite light curve, the median time between samples is 7.2 s. Finally, in Fig. [7] 
we show a time-binned version of the composite light curve to allow a visual comparison with the best previously-measured light 
curves, at optical and mid-infrared wavelengths. 

To determine the "allowed range" of each parameter — or, more precisely, the a posteriori joint probability distribution of all 
the parameter values — we employed the Markov chain Monte Carlo (MCMC) technique (see, e.g., Winn et al. 2007; Burke et 
al. 2007). We produced 8 chains of length 7.4 x 10 6 using a Gibbs sampler and a Metropolis-Hastings jump condition such 
that each parameter had an effective chain length of roughly 2 x 10 5 . This was accomplished by adjusting the scale of the 
jump-function distribution such that the probability of accepting a jump is roughly uniform and equal to approximately 40% 
across all parameters. To establish initial estimates of parameter uncertainties, a preliminary Monte Carlo bootstrap analysis was 
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FIG. 5. — NICMOS photometry (1.1-2.0 fim) of 4 transits of HD 149026b, after correcting for systematic effects. In each panel, the solid line shows the 
best-fitting model. 

performed, based on the Levenberg-Marquardt least squares minimization; then, the MCMC initial conditions were drawn from 
normal distributions with widths equal to five times these initial error estimates. The first 25% of each chain was trimmed, and 
then all the chains were concatenated. Each parameter had a Gelman & Rubin ( 1 992) R statistic smaller than 1 .0 1 , a sign of good 
convergence of the posterior parameter distributions. For each parameter, the values at each link of the chain were sorted. To 
describe the results, we report the median (50%) value, along with the interval between the 15.85% and 84. 15% levels (the 68.3% 
confidence interval). The results are given in Column 2 of Table ([T). 

3.1. Results from NICMOS photometric analysis 

For the orbital inclination, we find i = 84. 55^ ^ deg. This is in agreement with the independent estimate of i = 85.4^9 g deg by 
Nutzman et al. (2008), using the 8 /im channel of the Infrared Array Camera (IRAC) aboard the Spitzer Space Telescope. For the 

This is also in agreement with the results of Nutzman et al. (2008), who 
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FIG. 6. — NICMOS transit light curve (1.1-2.0 /im) of HD 149026b. The data from 4 transits have been superimposed, after correcting for systematic effects. 
The solid line shows the best-fitting model. 

found a/R+ = 6.20^ ??. The new result is more precise, which (as we will show in § |4j leads to tighter constraints on the stellar 
mass and radius. This is important because the uncertainties in the stellar properties have been the limiting factors in the analysis 
of this system. Based on the preceding results, we find the impact parameter (defined as b = acosi/R+) to be 0.57 l^^g. This 
is the tightest such constraint that has been achieved for HD 149026b. The earliest measurements of the impact parameter were 
consistent with zero (Sato et al. 2005, Charbonneau et al. 2006, Winn et al. 2008b), a situation leading to strong degeneracies 
among the transit parameters (Carter et al. 2008). More recently, Nutzman et al. (2008) found b = 0.62^ -^j, consistent with the 
new and more precise result. The increased precision in a //?* and b is a direct consequence of the improved precision with which 
the ingress (and egress) duration is known (Carter et al. 2008). In this sense, the greatest value of the NICMOS data is in the 
good coverage of the ingress and egress phases. 

The enhanced precision of the NICMOS data relative to previous data sets does not lead to a correspondingly enhanced 
precision in the planet-to-star radius ratio. This is because we allowed the time series from each orbit to have its own adjustable 
flux multiplier. Consequently, all of the information about R p /R+ comes only from those orbits that span an entire ingress or 
egress event. However, the value of R p /R+ that we derive is at least comparable in precision to previous determinations. We find 
Rp/R-k = 0.05416^01]™^. Interestingly this is larger by 2cr than the values derived previously, which were based on optical and 
mid-infrared data. Winn et al. (2008b) found R p /R+ = 0.049 1^;^ based on Stromgren (b+y)/2 photometry, and Nutzman et 
al. (2008) found R p /R* = 0.05158 ± 0.00077 based on 8 fxm photometry. 

Since we do not understand all of our noise sources with a physically-grounded model, we cannot rule out the possibility that 
the discrepancy between our result and the previous results is due to a faulty model of the systematic effects. The culprit would 
probably need to be the polynomial function of orbital phase. We do find that our result is unaffected if we replace the polynomial 
function of (f> with trigonometric functions, as mentioned previously; and we also find that our results are unchanged if we use a 
linear limb darkening law or fix the quadratic limb darkening parameters to those tabulated by Claret (2000). These tests do not 
prove that our results are valid but they do suggest that our analysis procedure is robust to changes in the functional form of the 
model. However, to the extent that the intra-orbital variations are not strictly repeatable within a given visit, our model would 
produce biased results. Fig.[H]shows the data after dividing by the optimized values of f x a (the orbit-specific flux multipliers) and 
dividing by the optimized transit model. The purpose of the divisions is to isolate the intra-orbital systematic effects, which do 
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FIG. 7. — Comparison of the best available transit light curves of HD 149026. Top panel: optical [Stromgren (b+y)/2] photometry from Sato et al. (2005) and 
Winn et al. (2008b), with a time sampling of 8.6 s and an rms residual of 2017 ppm. Middle panel: near-infrared [1.1-2.0 fim] photometry from this work, with 
a time sampling of 7.2 s and an rms residual of 440 ppm. Bottom panel: mid-infrared [8 fim] photometry from Nutzman et al. (2008), with a time sampling of 
7.4 s and an rms residual of 1854 ppm. The solid lines show the best-fitting model, which was calculated as described in §|5] 

appear consistent among the orbits within a given visit. 

Another possibility for the discrepancy in the transit depth between our near-infrared result and the previous optical and mid- 
infrared results is differential absorption due to constituents in the outer, optically-thin portion of the planet's atmosphere. This 
interpretation is the basis of the "transmission spectroscopy" technique for identifying constituents of exoplanetary atmospheres 
pioneered by Charbonneau et al. (2002). We consider this possibility at some length in sections [5] and [7] 

We find that the center-to-limb variation is less pronounced than was expected based on the tabulated limb darkening coeffi- 
cients of Claret (2000). Fig.|9]shows the confidence contours in the u\, U2 parameter space. The open square corresponds to the 
tabulated values for H band (for a star with T e ff = 6250 K, logg* = 4.5, [Fe/H] = 0.3, matching the properties of HD 149026). 
The tabulated values are excluded with >95% confidence. 

Two quantities intrinsic to the star and planet that may be calculated directly in terms of observables are the surface gravity of 
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FIG. 8 . — Isolation of the intra-orbital variations. The flux time series has been divided by the optimized flux multipliers f% and by the optimized transit model. 
The remaining variation appears to present a consistent pattern among all orbits within a given visit, as assumed in our model. The solid line is the optimized 
model. Each column shows data from all orbits of a given visit. Each row shows orbits arranged from first to last in rows from top to bottom, respectively. 

the planet, and the mean density of the star. The surface gravity of the planet is calculated as (Southworth et al. 2007) 



2tt 



K 



g P -- 



(9) 



P {R p / a) 2 smV 

where K is the semiamplitude of the stellar radial-velocity variation (43.3 ± 1.2 m s" 1 ; Sato et al. 2005), and R p /a and sin / are 
derived from our MCMC analysis. We find logg^ = 3.132+QQ3J where g p is in cgs units. The mean stellar density p+ is calculated 
as (Seager & Mallen-Ornelas 2003, Sozzetti et al. 2007) 



3ir ( a 

gp^KpZ 



-Pp 



(10) 



where p p is the mean density of the exoplanet. We may neglect the correction term involving the planetary density as p + ~ p p , 
(Rp/R*) 3 ^1 x 10" 4 and the fractional uncertainty in a/R± is larger than the fractional uncertainty in Rp/R* by a factor of 2. 
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FIG. 9. — Results for the limb-darkening parameters u\ and ui- The contours are the 68% and 95% confidence regions as determined by the MCMC analysis 
of the photometric time series. The solid square is the minimum-x 2 solution, u\ = 0.0, 112 = 0.1789. The open square marks the tabulated values of Claret (2000) 
for H band (h, = 0.044, u 2 = 0.344). 

4. STELLAR PARAMETERS 

The basic inputs to models of the planetary interior are the planetary mass M p and radius R p , in units of grams and kilometers, 
or in units of Jupiter's mass and radius. Transit photometry and Doppler velocimetry alone do not determine these quantities. 
Additional information about the star must be introduced. Several techniques for estimating the stellar mass M 4 and radius 
were reviewed by Winn et al. (2008b). We chose to estimate and using stellar-evolution models that are constrained by the 
best available, relevant, observable properties of the star: the mean density 0.497^' JJ^ g cm" 3 determined from our light-curve 
analysis, the absolute magnitude My = 3.65 ±0.12 derived from the Hipparcos parallax and apparent magnitude [tt= 12.59 ±0.70 
mas, V = 8.15 ±0.02; van Leeuwen (2007)], effective temperature [T e g = 6160 ±50 K, a weighted mean of the results from Sato 
et al. (2005) and Masana et al. (2006)], and metallicity [0.36 ±0.08, from Sato et al. (2005) with a more conservative error bar]. 
We chose not to use the spectroscopically-determined stellar surface gravity [logg* = 4.26 ±0.07; Sato et al. (2005)] because 
the photometrically-determined value of p+ provides an effectively tighter constraint, and because spectroscopically-deter mined 
surface gravities have been found to be susceptible to systematic error (see, e.g., Winn et al. 2008a. IJohns-Krull et al.l d2008l) ). 

Following the procedure of Torres, Winn, & Holman (2008), we employed Yonsei-Yale stellar models 6 (Yi et al. 2001, Demar- 
que et al. 2004). Model isochrones were interpolated in both age and metallicity, for metallicities [Fe/H] ranging from 0.28 to 
0.43 and for ages ranging from 0.1 to 14 Gyr, in steps of 0.1 Gyr. Fig. [TOl shows several of these theoretical isochrones, along 
with some of the observational constraints. The upper left panel illustrates the constraint due to the spectroscopically-determined 
surface gravity, even though we did not actually apply that constraint, as explained above. It is evident that the constraint due to 
p+ is stronger. 

The isochrones were interpolated to provide a fine grid in stellar mass (with a step size of 0.005 M Q ). We then assumed that 
the likelihood of each point on the interpolated isochrones was proportional to exp(-\l/2), where 

<2 _ / A[Fe/H] \ 2 | /Ar eff \ 2 | / AM V \ 2 | / A Pir 

V CT[Fe/H] / V °Tm J V a Mv J V °>. 

and the A quantities denote the differences between the observed and calculated values. The asymmetry in the error distribution 
for p± was taken into account. Additionally, the likelihood was taken to be proportional to an Salpeter initial mass function, 

6 We chose the Y 2 models mainly because of the convenient form in which they are publicly available. Other stellar-evolutionary models are available, and 
other investigators have examined the sensitivity of results such as ours to the choice of model. For HD 149026 in particular, Southworth (2008) found that the 
Y 2 models and independent models by Claret (2007) gave results for and that agreed to within 1%. Another set of publicly available models, the Padova 
models of Girardi et al. (2000), are not computed for the high metallicity observed for HD 149026; but at zero metallicity, at least, both Torres et al. (2008) and 
Southworth (2008) found that the Padova models give results for the stellar mass and radius that are also within 1% of the Y 2 results. 
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Fig. 10. — Stellar-evolutionary model isochrones, from the Yonsei-Yale series by Yi et al. (2001). The points and shaded boxes represent the observationally- 
determined values and la errors. Here, surface gravity is determined spectroscopically (Sato et al. 2005), My is derived from Hipparcos parallax and V 
magnitudes, and is determined photometrically from the transit light curve. Isochrones are shown for ages of 1 to 13 Gyr (from left to right) in steps of 1 Gyr 
for a fixed stellar metallicity of [Fe/H]= 0.36. 

£(M) oc M" (1+l) with x = 1 .35 (Salpeter 1955). The joint probability function, P, was taken to be proportional to the likelihood, 
viz., 

P(^,M + ,r eff ,log^,M v ,[Fe/H],p + ,Age) cxexp(- X 2 j2). (12) 
For a given parameter Xq from this list, we calculated the cumulative distribution function (CDF) by numerically evaluating 

/X poo />OC 

dX / dX 1 - - dX N P{X ,X X , ■ ■ ■ ,X N ). (13) 
oo J— oo J —oo 

For each parameter, we record the values of x for which the CDF takes the values 15.85%, 50%, and 84. 15%. The 50% level (the 
median) is reported as the "best-fit value" and the interval between the 84.15% and 15.85% levels is reported as the 68.3% (lex) 
confidence interval. 7 

Table (Q]) reports the best-fit stellar parameters and confidence intervals. We find the stellar radius to be = 1.541^ j]^ ^o- 
This is larger than (but in agreement with) the previous estimates of = 1.46 ±0.10 Rq by Sato et al. (2005), and R+ = 
1 .497 ± 0.069 Rq) by Nutzman et al. (2008). By combining the derived distribution for R+ with the photometrically-determined 
distribution for Rp/R^, we find the planetary radius to be R p = 0.813±q q25 ^Jup- This is larger than any previous result. Using 
optical photometry, Sato et al. (2005) found 0.725 ± 0.050 R Jup , Charbonneau et al. (2006) found 0.726 ± 0.064 R ]up , and Winn 
et al. (2008b) found 0.71 ± 0.05 R ]up . Using mid-infrared photometry, Nutzman et al. (2008) found 0.755 ± 0.040 R ]up . It is 
important to note that these determinations were not wholly independent, and therefore should not be combined into a weighted 
average. They all used many common inputs for the stellar properties, and the analyses of optical photometry all included a 
common subset of at least 3 light curves. 

7 Although our procedure was inspired by the work of Torres et al. (2008) and is similar in almost all respects, there is one significant difference. The best 
fit values reported by Torres et al. (2008) were those that minimized \~ as m our analysis. The difference is that Torres et al. (2008) estimated the la errors in 
the stellar properties based on the total span of the calculated values that gave agreement within la with the observables. Effectively, they assumed a uniform 
error distribution for each observable, rather than a Gaussian error distribution as we have done. Consequently, our method produces smaller error intervals in 
the stellar properties. Caution would dictate that larger error intervals are desirable, especially since we are relying on the theoretical isochrones that surely have 
some unaccounted-for systematic errors. However, using a uniform error distribution for the observables is an arbitrary way to inflate the output errors, and the 
true error distribution for the observables is probably closer to Gaussian. For these reasons we chose our approach and emphasize the caveat that our results place 
complete trust in the Y" isochrones. 
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As mentioned previously, we did not apply any constraint to the models based on the spectroscopically determined value of 
logg*. However, given our results for and R+ we computed the implied value of logg*, finding logg* = 4.189^ qjj. This in 
agreement with, and is more precise than, the spectroscopically-determined value of logg* = 4.26 ± 0.07 (Sato et al. 2005). 

5. JOINT ANALYSIS WITH OPTICAL AND MID-INFRARED LIGHT CURVES 

Transit observations of HD 149026b have now been made at optical wavelengths (Sato et al. 2005, Charbonneau et al. 2006, 
Winn et al. 2008b), near-infrared wavelengths (this work), and mid-infrared wavelengths (Nutzman et al. 2008). In this section 
we repeat our analysis on all of these data, in order to bring all of these data to bear on the determination of the system parameters, 
while seeking possible wavelength variations in the planet-to-star radius ratio. 

It is reasonable to require consistency across these data in the parameters relating to the orbital configuration of the transit, 
such as the inclination angle and normalized semi-major axis. However, the inferred planet-to-star radius ratio is a wavelength- 
dependent quantity, depending on the opacity of the exoplanetary atmosphere and the emergent flux from the planetary nightside 
(which is expected to be unimportant). With this in mind, we performed a joint analysis of all of the data, requiring consistency 
in cos/ and a/R+ but allowing R p /R* to take separate values for each of the three types of data: optical, near-infrared, and 
mid-infrared. 

Specifically, we fitted our NICMOS data, the 8 /im IRAC time series of Nutzman et al. (2008), and the 8 light curves obtained 
in the Stromgren (b + y)/2 band by Sato et al. (2005) and Winn et al. (2008b). Our photometric model for the NICMOS data, 
including the associated systematic effects, has already been described. For the data sets presented by other authors, we followed 
those authors' prescriptions to account for systematic errors. For the (b+y)/2 data, we corrected the data by allowing the out- 
of-transit flux to be a linear function of time. For the 8 /im data, we modeled the time-variable sensitivity of the detector (the 
"ramp") as a multiplicative correction, / sys = ao + «ilog(f-fo) + «2log 2 (/-fo), where to is the time immediately prior to the start 
of the observation. 

We performed an MCMC analysis of this joint data set. The free parameters relating to the transit model were the three values 
of R p /R+ (corresponding to the ratios measured at approximately 0.5 [im, 1.5 [im, and 8.0 /im); the geometric parameters cos; 
and a/R+; the quadratic limb-darkening coefficients for the NICMOS light curve; the linear limb-darkening coefficients for the 
optical and infrared light curves (for which the precision of the data do not justify the more accurate quadratic law); and the 
mid-transit times for the NICMOS data and the IRAC data. 8 We also fitted for the ramp-correction terms for the IRAC data and 
the parameters relating to the flux offsets and intra-orbital variations for the NICMOS data. Six chains of length 9 x 10 6 were 
created, representing approximately 2 x 10 5 correlation-lengths per parameter. These were concatenated after removing the first 
25% of each chain. The Gelman-Rubin R statistic was smaller than 1.01 for each parameter. We then repeated the analyses that 
were described in §[3]and §|4]to determine the stellar, planetary, and orbital parameters, based on this joint analysis. The results 
are tabulated in Table (Q~|). 

The results for the geometric parameters are hardly changed from the NICMOS-only analysis, a reflection of the greater 
precision of the NICMOS light curve. The planet-to-star radius ratio was found to be larger for the NICMOS data than for the 
other bandpasses, as was already evident from the comparison of our NICMOS-only analysis to previously published analyses. 
Fig- H3 shows the variation in (R p /R+) 2 with wavelength. The quantity (R p /R+) 2 is essentially the transit depth, or fractional 
loss of light during the total phase of the transit, after "removing" the effects of limb darkening. The radius ratios that we 
derive for the mid-infrared and optical data are in agreement with those reported previously. For the IRAC data we find R p /R* = 
0.05188^:8oo86 as compared to the value 0.05158 ± 0.00077 reported by Nutzman et al. (2008). For the (b+y)/2 data we find 
R p /R+ = 0.05070;;S;[J[|JJ| as compared to the value 0.0491+g ^ | found by Winn et al. (2008b). The precision in the optical R p /R* 
has been increased because the NICMOS data pins down all of the other parameters that are correlated with R p /R+. 

6. EPHEMERIS AND TRANSIT TIMING 

The NICMOS-only analysis resulted in the measurement of four distinct mid-transit times, with uncertainties smaller than 
45 s. These are given in Table Q. We pooled together all of the independent measurements of mid-transit times from Winn et 
al. (2008b), Nutzman et al. (2008), and this work, to derive a new transit ephemeris. We fitted the times to a linear function of 
the integral epoch E, 

T C (E) = T C {Q) + EP (14) 

where P is the period and T c (0) is the mid-transit time at some fiducial epoch. The results were T c (0) = 2454456.78761 ±0.00014 
HJD and P = 2.875891 1 ± 0.0000025 days. The linear fit had x 2 = 20.16 and 14 degrees of freedom. This is a marginally 
unacceptable fit. The formal probability to find a value of \ 2 this large is 15%. Further transit observations are needed to 
distinguish the possibilities of a genuine period variation, a statistical fluke, and underestimated timing errors. Fig.QT]shows the 
O-C (observed minus calculated) timing diagram. 

7. DISCUSSION 

We have presented observations of four transits of HD 149026b at near-infrared wavelengths with the HST NICMOS detector. 
The NICMOS data place the strongest constraints yet on the geometrical system parameters. In particular, the increased precision 
of the measurement of the normalized semi-major axis (a/R+) leads to an improved estimate of the mean stellar density, which 
was then coupled with stellar evolution models to constrain the stellar mass and radius. Improved knowledge of the stellar mass 

8 To keep the number of parameters as small as possible, we did not vary the optical mid-transit times or baseline correction parameters at this stage, having 
found that they are uncorrelated with the other parameters of interest. 
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FIG. 1 1 . — Transit-timing variations for HD 149026b. The differences between observed and calculated mid-transit times are plotted. The last 4 points represent 
the new NICMOS observations. The dashed lines show the Icr range in the calculated times according to the linear ephemeris presented in §|6] 

and radius leads to greater precision in the planetary mass and radius. We have found a larger stellar radius, and a larger planet- 
to-star radius ratio, than previous estimates. As a result of these two factors, we have also found the planetary radius to be larger 
than previously thought. The planet has "grown" by about 7%. 

Despite this increase, our results are still consistent with the contention that HD 149026b is highly enriched in heavy elements. 
It is still smaller than expected for a hydrogen-helium planet with the given mass and degree of stellar irradiation (Burrows et 
al. 2007). For comparison, the tabulated models by Fortney et al. (2007) predict a 1.3 Rj np hydrogen-helium HD 149026b at an 
age of 1 Gyr. A variety of models have been developed to estimate the heavy-element content of HD 149026b (Sato et al. 2005, 
Fortney et al. 2006, Ikoma et al. 2006, Burrows et al. 2007), most of which suppose that the metals are confined to an inner 
core of material beneath a hydrogen-helium envelope. Other physical considerations in these models include the equation of 
state for heavy elements at core pressures, atmospheric opacities and the upper boundary condition where energy is delivered 
from the star. To determine a revised estimate for the heavy-element content, we used the tabulated models provided by Sato et 
al. (2005) and Fortney et al. (2007), and interpolated the tabulated results as appropriate for the planetary radius, planetary mass, 
and degree of irradiation that follow from the parameters determined from the NICMOS data. We find a core mass in the range 
of 45-70 M®, depending on assumed stellar age and core density. Thus, the interpretation of the planet as highly enriched is 
unaffected, although the required amount of enrichment is slightly reduced. 

It is also interesting that the planet-to-star area ratio, (R p /R+) 2 , was found to be 2-3 er larger in the NICMOS band (1.1-2.0 /im) 
than in the optical band (0.45-0.55 /im) or mid-infrared bands (6.5-9.5 /im), while the results for the latter two bands are in 
agreement. Caution dictates that this discrepancy should not be over-interpreted. It is possible that the discrepancy is at least 
partly the result of unresolved systematic errors in any of the data sets. We have already noted that the noise in the NICMOS data 
exceeds the photon noise level by a factor of 2, and is not well understood. 

However, it is also worth considering that this wavelength-dependent variation represents selective absorption by constituents 
in the outermost layer of the planet's atmosphere. Molecules with strong absorption bands at near-infrared wavelengths would 
cause the transit to appear deeper at those wavelengths. Strong bands are expected for the common molecules CO, H2O, and 
CH 4 (Brown 2001, Hubbard et al. 2001, Seager & Sasselov 2000). If this were the case, then a detailed analysis of the NICMOS 
spectrophotometry — breaking it down into smaller wavelength bins, as opposed to summing the entire first-order spectrum — 
might be used to identify some constituents of the planet's atmosphere. In addition, more care would be needed in choosing 
which radius to use in the comparison with models of the planet's interior. It is beyond the scope of this paper to analyze the 
wavelength dependence of the transit depth across the NICMOS band, or to compute a realistic atmospheric model to see if 
the contrast between the optical, near-infrared, and mid-infrared results can be accommodated. We can, however, perform an 
order-of-magnitude calculation to check on the plausibility of this interpretation. 

Let z(t) be the height in the planet's atmosphere at which the optical depth is t for a path from the star to the observer, as 
diagrammed in Fig. [T2] This height is measured relative to an atmospheric base radius Rq, where the planet is optically thick 
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at all relevant wavelengths. The height z(t) depends, in part, on the wavelength-dependent opacity and the density profile of 
the atmosphere. We define R P (X) as Rq + z(t = 1) and S = (R p /R+) 2 . If we assume that z(l) <C ^o, then the transit depth is 
approximately linear in z(l): 
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Next, we consider the difference in S as measured in two distinct wavelength bands: 
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where we have defined the height difference Sz = Zi(l)~Z2(l)- The height difference Sz reflects differing levels of absorption in 
the two bands. Solving for Sz in Eqn. ( [Tol l we find 
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FIG. 12. — Illustration of wavelength-dependent absorption. Shown are some rays that skirt the planetary atmosphere on their way to Earth. At some height 
z(t) above the fiducial radius Rq, the ray (parameterized by path length s) has an optical depth of r (solid line). Light that follows paths with z < z(j = 1) are 
mainly absorbed. The height z(r = 1) corresponding to optical depth of unity depends on wavelength, giving rise to a wavelength-dependent transit radius. 

For HD 149026b, to evaluate the idea that the larger near-infrared measurement of <5 is due to molecular absorption, we assume 
that Rq is the optically-derived radius (Rq = 0.757 R] up ) because the optical spectrum is expected to show comparatively weak 
absorption features (Brown 2001, Seager & Sasselov 2000). Using Eqn. dPTl ). the 0.035% difference in S between near-infrared 
and mid-infrared wavelengths implies Sz ~ 2500 km (4% of Rq). To judge if this is realistic, we are interested in expressing Sz 
in units of the pressure scale height, for which an order-of-magnitude expression is H = kT / [im p g p , where T is a representative 
atmospheric temperature, /i is the mean molecular weight of the atmosphere, m p is the proton mass and g p is the surface gravity. 
Using the surface gravity that was determined from our analysis of the optical light curve (g p = 1535 cm s" 2 ), and assuming a 
H2-He atmosphere with T = 2300 K [as measured at 8 /im by Harrington et al. (2007)], we find H k, 530 km, and Nh ~ 5. If 
instead we use the planet's predicted temperature at thermal equilibrium with the incident stellar radiation (T = 1700 K), we find 
H w 400 km and Nh ~ 6. 

If we assume further that the absorbers have an exponential density profile, 



p(z) = p(0)exp(-|) , 



(18) 



where h is the density scale height, then we may express z(r) in terms of the opacity a of the absorbing molecules, as follows. By 
integrating the optical depth t across the optical path at height z^Rq (as illustrated in Fig.[T2l. where a change of path length 
ds results in a change in optical depth dr as 
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In general, h may be different for each atmospheric constituent. If we assume that the components are uniformly mixed throughout 
the atmosphere, then h is independent of composition, and the difference between two heights Zi(r) and Z2(r) at two distinct 
wavelength bands with different opacities a\ and 02 can be written independently of the optical depth, as 



Sz = zi -Z2 ~h\n [ — 

, a 2 



(21) 
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If we assume further that the temperature scale height is large compared to the pressure scale height, then h w H and Nh ~ 
ln(oi/<7 2 ). 

We may now judge the plausibility of this interpretation with reference to the typical opacities and widths of molecular absorp- 
tion features. For strong molecular bands and atomic lines, the ratio of the in-band opacity to the nearby continuum opacity may 
be as large as 10 4 (Brown 2001, Seager & Sasselov 2000), yielding a maximum height difference of Nh ~ 10 scale heights within 
the absorption band. The NICMOS band from 1.1-2.0 /im includes strong rotation-vibration molecular absorption bands due to 
H2O, CO and CH4. An example of a very strong absorption band is a water band centered at 1.4 /im, spanning approximately 
10% of the effective filter width. If we assume that this is the dominant spectral feature in this band, then the result of Nh ~ 5 
across the entire bandpass translates into Nh ~ 50 within the bandpass of the absorption feature. This is larger than the criterion 
Nh ~ 10 mentioned above. 

Therefore this interpretation seems to require significantly more opaque or broader-band absorption features than are seen in 
the models. In one sense the result of the order-of-magntiude calculation is discouraging, as it may make it seem more likely that 
the discrepancy in depths is due to systematic errors. On the other hand, if the noise were well-understood and the discrepancy 
could be confidently proclaimed, then it would be the sign of new and interesting atmospheric physics that is not described in the 
standard models. Some priorities for progress on this issue include an examination of the wavelength-dependence of the transit 
across the NICMOS band, and the observation of the system with other NICMOS grisms, which are reputed to be more stable 
than the G141 grism used here. 
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TABLE 1 

System Parameters of HD 149026. 



Epoch Mid-transit time [HJD] Error 

2454456.78751 0.00030 

1 2454459.66379 0.00023 
17 2454505.67688 0.00034 
31 2454545.94133 0.00047 



TABLE 2 

Mid-transit times, based on the NICMOS data. 
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FIG. 13. — The planet-to-star area ratio, (R p /R^) 2 , as a function of observing wavelength, based on a joint fit to the NICMOS data, the (b+y)/2 data of Winn 
et al. (2008), and the IRAC data of Nutzman et al. (2008). The horizontal error bars show the approximate wavelength range of each bandpass. 



